library(ggplot2)
require(data.table)
library(ncdf4)

#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\OneDrive - NORCE\\NFood_replication\\"

crops = c("Corn", "Wheat", "Soybean", "Cotton", "Rice", "Sugarcane", "Tropical Corn", "Tropical Soybean")
case = c("RCP45","RCP45_85LU","RCP85","RCP85_45CO2","RCP85_SAI","RCP85_MSB","RCP85_CCT")


#==================Final: plot original and soy-corrected co2 response in one plot
#see: script "Data.ED.Fig.3.R"
load(paste0(wd1,"/data/RData/Fig.ED3.CO2_effect.RData"))

library(forcats)
#jpeg(paste0(getwd(),"/crop_figures/final/Fig.Ext4.co2_only_effects_pct_quadratic.jpeg"), units = "in", width = 6, height = 4, res = 600)
p <- ggplot(df.co2_cor, aes(co2.dif, effect, group=Crop)) + #reorder the group sequence based on median
  geom_line(aes(linetype = fct_reorder(Crop, effect, .fun=median, .desc = TRUE), #color = Crop, 
                colour = fct_reorder(Crop, effect, .fun=median, .desc = TRUE))) + #must reorder linetype, colour, shape consistently
  geom_point(aes(co2.dif,yield.dif, group=Crop, shape = fct_reorder(Crop, effect, .fun=median, .desc = TRUE),
                 colour = fct_reorder(Crop, effect, .fun=median, .desc = TRUE) )) +
  scale_linetype_manual(values=c(1,1,2,4,3,6,5))+ 
  scale_shape_manual(values=c(16,16,1,2,17,8,0)) +
  scale_colour_manual(values=c("#000000","#D55E00","#000000","#000000","#000000","#000000","#000000"))+ #,  guide = FALSE)+
  scale_x_continuous(limits=c(0, 400), breaks =c(0, 100, 200, 300, 400)) 
  # guides(colour=guide_legend(title=NULL, nrow=2, byrow=T,  keywidth = 2, keyheight=1, label.position = "right", label.hjust = 0,
  #                            override.aes = list(linetype=1, shape=c(16,16,1,2,17,8,0), alpha=0 )))+
  # shape =guide_legend(override.aes = list(size = 10, shape=c(19,19,19,19,4), alpha=0.6))) +
##remove grid and background or remove all
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
p <- p + xlab(expression('CO'[2]~'change (ppm)')) + ylab("Yield change (%)") #use * instead of paste0 in expression
p <- p + theme(legend.title = element_blank(), legend.text = element_text(size=10)) + 
  theme(axis.text = element_text(size=9), axis.title=element_text(size=10))

print(p)
#dev.off()
#ggsave(paste0(wd1,"/figures/extended/Fig.Ext3.svg"), device ="svg", units = "in", width = 6, height = 4, dpi = 600)
ggsave(paste0(wd1,"/figures/extended/ED_Fig3.jpg"), device ="jpg", units = "in", width = 6, height = 4, dpi = 350)













